import math

import networkx

from Interface import Calculator


def cm(n, x, N, M, e=1e-6):
    '''cumulative hypergeometric distribution
    Given N balls, M of which are black and the rest are white, 
    what is the probability C(n, x,N,M) that out of n balls
    selected uniformly at random without replacement, at most x are black
    e is error
    result should less than 1
    '''
    if x == n:
        return 1.0
    
    def inv_jm(n, x, N, m):
        return (1-x/(m+1))/(1-(n-1-x)/(N-1-m))

    n = float(n)
    x = float(x)
    N = float(N)
    M = float(M)
    s = 1.0

    for k in xrange(int(x), int(M)-1):
        s = s * inv_jm(n, x, N, k) + 1.0

    ak = s
    bk = s
    k = M-2
    epsk = 2*e

    while k< (N-(n-x)-1) and epsk > e:
        ck = ak/bk
        k += 1
        jjm = inv_jm(n, x, N, k)
        bk = bk * jjm + 1.0
        ak = ak * jjm
        epsk = (N-(n-x)-1-k)*(ck-ak/bk)

    return 1-(ak/bk-epsk/2)


class PPCalculator(Calculator):
    
    def __init__(self, db, **keywds):
        Calculator.__init__(self, db)
        self.raw_result = networkx.Graph()
        
    def settings(self):
        pass
    
    def calculate(self, taxon=None):
        taxon_set = self.db.get_taxons()
        N = len(taxon_set)
        if taxon == None: # get a random taxon
            taxon = taxon_set.pop()
        profile = self.db.get_profiles(taxon)
        id_list = profile.keys()
        
        for i in xrange(0, len(id_list)-1):
            id1 = id_list[i]    
            profile1 = profile[id1]
            m = len(profile1)            
            if m == 0: continue
    
            for j in xrange(i+1, len(id_list)):
                id2 = id_list[j]
                profile2 = profile[id2]
                n = len(profile2)
                if n == 0: continue
    
                k = len(profile1.intersection(profile2))
                if k == 0: continue
            
                p = 1-cm(n, k, N, m)
                self.raw_result.add_nodes_from([id1, id2])
                self.raw_result.add_edge(id1, id2, p)
    
    def write(self, handler):
        # do nothing
        handler.close()
    

def _test():
    print cm(3,3,100,70)
        
        
if __name__ == '__main__':
    _test()
    
        
        
#END
    
        